R Notebook

library(readr)
library(ggplot2)
library(dplyr)
library(rms)
library(synthpop)
source('./helper functions.R')

Reading in original dataset, performing multiple imputation and writing out imputed dataset. This chunk commented out because we’ve already done this and saved the resulting file as our starting point.

Read in dataset from Brinati paper. Dataset already been imputed based on code above.

require(readr)
brinati_covid_study_v2.imputed <- read_csv("data/brinati-covid_study_v2_imputed.csv",
  col_types = cols(GENDER = col_factor(levels = c("M",
  "F")), SWAB = col_factor(levels = c("0",
  "1"))))

p<-plot_model(glm.orig.fit)
p<-plot_model(glm.orig.fit)
p<-ApplyFigureThemeLargeFontOnly(p + ylim(.1, 2.5) ) 

Scale for 'y' is already present. Adding another scale for 'y', which will replace the existing scale.

Draw ROC curve

p <- ff %>% mutate(fpr=1-Specificity) %>%
  ggplot(., aes(x=fpr, y=Sensitivity)) + geom_point(aes(frame=Cutoff)) +
    xlab('False Positive Rate (1-Specificity)') + ylab('True Positive Rate\n(Sensitivity)')

Ignoring unknown aesthetics: frame

ggplotly(p)

Draw calibration

Calculate hosmer lemeshow statistics

require(performance)
hl.org <- performance_hosmer(glm.orig.fit)
hl.org

# Hosmer-Lemeshow Goodness-of-Fit Test

  Chi-squared: 5.333
           df: 8    
      p-value: 0.722

Summary: model seems to fit well.

Create smoothed cal curve from Van hoorde et al

p<-CreateSmoothedCalCurvePlot(predicted, observed)
p<-ApplyFigureThemeLargeFontOnly(p)
SaveStdSquareFigure(p, 'figs/smoothed-cal-curve-orig-data.png')
p

Boot strap performance in the original dataset

Generate new a synthetic dataset with a prevalence of 10%

#Assess performance of the original model on new prevalence 10 data
assessPerf(predicted = plogis(predict(glm.orig.fit, newdata=as.data.frame(brinati.syn.prev10))), 
           observed = brinati.syn.prev10$SWAB)

          Dxy       C (ROC)            R2             D      D:Chi-sq           D:p             U      U:Chi-sq 
 6.462721e-01  8.231360e-01 -1.737270e+00 -6.089445e-01 -1.688955e+02            NA  7.433448e-01  2.093932e+02 
          U:p             Q         Brier     Intercept         Slope          Emax           E90          Eavg 
 0.000000e+00 -1.352289e+00  2.037312e-01 -2.471404e+00  6.986839e-01  5.722161e-01  5.415682e-01  3.070988e-01 
          S:z           S:p 
 3.691032e+00  2.233457e-04 

Generate new a synthetic dataset with a prevalence of 63% – close to the original dataset

#Assess performance of the original model on new prevalence 10 data
assessPerf(predicted = plogis(predict(glm.orig.fit, newdata=as.data.frame(brinati.syn.prev10))), 
           observed = brinati.syn.prev10$SWAB)

          Dxy       C (ROC)            R2             D      D:Chi-sq           D:p             U      U:Chi-sq 
 7.079320e-01  8.539660e-01 -7.403297e-01 -3.535701e-01 -9.764606e+01            NA  5.720065e-01  1.615898e+02 
          U:p             Q         Brier     Intercept         Slope          Emax           E90          Eavg 
 0.000000e+00 -9.255766e-01  1.923112e-01 -2.097882e+00  7.796140e-01  4.378514e-01  4.367189e-01  2.627114e-01 
          S:z           S:p 
 3.977260e+00  6.971398e-05 

Boot strap performance of a complete re-moel of a prevalence 10% dataset

2. Look at performance as a function of prevalence

prevalences = seq(0.01, 0.99, 0.01)
dfs_prevs <- lapply(prevalences, function (prev) brinati.syn.factory(prev = prev))
dfs.prevs.perf <- lapply(dfs_prevs, function(df) {
  assessPerf(predicted = plogis(predict(glm.orig.fit, newdata=as.data.frame(df))), observed = df$SWAB)
})

Calculate the different validation measures

dfs.prevs.perf.df <- data.frame(matrix(unlist(dfs.prevs.perf), nrow = length(prevalences), byrow = TRUE))
colnames(dfs.prevs.perf.df) <- names(dfs.prevs.perf[[1]])
dfs.prevs.perf.df$Prevalence = prevalences
dfs.prevs.perf.df <- gather(dfs.prevs.perf.df, Measure, Value, Dxy:`S:p`, factor_key = TRUE)

Plot the distribution of validation measures

measures = c('C (ROC)', 'Brier', 'Intercept', 'Slope')
ggplot(dfs.prevs.perf.df %>% filter(Measure %in% measures), aes(x=Measure, y=Value,  color=Measure)) + geom_boxplot() +
  geom_point(data=data.frame(Measure = measures, Value = c(1, 0, 0, 1)), size=3, color='blue') + theme_bw()

ggplot(dfs.prevs.perf.df %>% filter(Measure %in% measures), aes(x=Prevalence, y=Value,  color=Measure, group=Measure)) + geom_line() + theme_bw()